In this document, we will fit all models regarding this project. The detailed preregistered analysis plan can be found in our OSF’s project.
Let’s load libraries and data.
# Ensures the package "pacman" is installed
if (!require("pacman")) install.packages("pacman")
# Install/Load packages
pacman::p_load("rio",
"here",
"metafor",
"dplyr",
"bayesmeta",
"brms")
remotes::install_github("stan-dev/cmdstanr")
# Ensure that packages are installed according to versions in renv's lockfile
renv::restore()
# Import data
d = rio::import(here::here("data", "data.xlsx"))
Let’s calculate the crude odds ratio based on the extracted data.
# Calculate crude log odds ratios
d_logOR =
metafor::escalc(
measure = "OR", # log odds ratio,
# Tocilizumab
ai = trt_events,
n1i = trt_total,
# Control
ci = control_events,
n2i = control_total,
data = d
) %>%
as_tibble() %>%
# Calculate standard error
dplyr::mutate(sei = sqrt(vi),
# Set order to use "tocilizumab" as the Intercept in brms
treatment = factor(treatment,
levels = c("tocilizumab",
"sarilumab",
"baricitinib")))
# saveRDS(d_logOR,
# here("output", "data", "effect_sizes.Rds"))
First, we will fit a single meta-analysis model per drug (tocilizumab, sarilumab, or baricitinib vs. control).
All Bayesian meta-analysis random-effect models are the same, and is described as:
\[ \begin{align*} y_i & \sim Normal(\theta_i, \sigma_i^2) \tag{Likelihood} \\ \theta_i & \sim Normal(\mu, \tau^2)\\ \\ \mu & \sim \operatorname{Normal}(0, 0.75^2) \tag{Priors} \\ \tau & \sim \operatorname{Log-Normal}(-1.975, 0.67^2) \\ \end{align*} \]
where \(y_i\) is the observed mean log odds ratio in study \(i\) for either tocilizumab, sarilumab, or baricitinib versus control treatment. We assume these effect sizes are normally distributed around the study-specific mean \(\theta_i\) along with a known sampling variance, represented by the observed \(\sigma_i^2\). We also assume \(\theta_i\) is drawn from a normal distribution where \(\mu\) is the average effect and \(\tau^2\) is the between-study heterogeneity.
## Formula
mf =
# https://bookdown.org/content/4857/horoscopes-insights.html#consider-using-the-0-intercept-syntax
formula(yi | se(sei) ~ 0 + Intercept + (1 | study))
## Priors
# Tau
informative = bayesmeta::TurnerEtAlPrior("all-cause mortality",
"pharmacological",
"placebo / control")
logmean = informative$parameters["tau", "mu"]
logsd = informative$parameters["tau", "sigma"]
# logmean
# -1.975
# logsd
# 0.67
priors_ma =
brms::prior(normal(0, 0.75), class = "b", coef = "Intercept") +
brms::prior(lognormal(-1.975, 0.67), class = "sd")
## Models
# Tocilizumab
ma_toci =
brms::brm(
data = d_logOR |> dplyr::filter(treatment == "tocilizumab"),
family = gaussian,
formula = mf,
prior = priors_ma,
sample_prior = TRUE,
control = list(adapt_delta = .95),
backend = "cmdstanr", # faster
cores = parallel::detectCores(),
chains = 4,
warmup = 2000,
iter = 4000,
seed = 123,
file = here::here("output", "fits", "ma_toci.Rds"),
file_refit = "on_change"
)
# Sarilumab
ma_sari =
brms::brm(
data = d_logOR |> dplyr::filter(treatment == "sarilumab"),
family = gaussian,
formula = mf,
prior = priors_ma,
sample_prior = TRUE,
control = list(adapt_delta = .95),
backend = "cmdstanr", # faster
cores = parallel::detectCores(),
chains = 4,
warmup = 2000,
iter = 4000,
seed = 123,
file = here::here("output", "fits", "ma_sari.Rds"),
file_refit = "on_change"
)
# Baricitinib
ma_bari =
brms::brm(
data = d_logOR |> dplyr::filter(treatment == "baricitinib", study != "RECOVERY Bari"),
family = gaussian,
formula = mf,
prior = priors_ma,
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
cores = parallel::detectCores(),
chains = 4,
warmup = 2000,
iter = 4000,
seed = 123,
file = here::here("output", "fits", "ma_bari.Rds"),
file_refit = "on_change"
)
ma_bari_sens =
brms::brm(
data = d_logOR |> dplyr::filter(treatment == "baricitinib", study != "RECOVERY Bari (No Toci)"),
family = gaussian,
formula = mf,
prior = priors_ma,
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
cores = parallel::detectCores(),
chains = 4,
warmup = 2000,
iter = 4000,
seed = 123,
file = here::here("output", "fits", "ma_bari_sens.Rds"),
file_refit = "on_change"
)
Now, we will fit Bayesian meta-regression models (one for tocilizumab vs. sarilumab data, another for tocilizumab vs. baricitinib).
The model can be described as:
\[ \begin{align*} y_i & \sim Normal(\theta_i, \sigma_i^2) \tag{Likelihood}\\ \theta_i & \sim Normal(\mu, \tau^2)\\ \mu &= \beta_0 + \beta_1 x_i\\ \\ \end{align*} \]
where \(y_i\) is the observed mean log odds ratio in study \(i\) for either tocilizumab, sarilumab, or baricitinib versus control treatment. We assume these effect sizes are normally distributed around the study-specific mean \(\theta_i\) along with a known sampling variance, represented by the observed \(\sigma_i^2\). We also assume \(\theta_i\) is drawn from a normal distribution where \(\mu\) is the average effect and \(\tau^2\) is the between-study heterogeneity. We will estimate \(\mu\) as a function of the estimated population effect size \(\beta_0\) and the moderator \(x\) multiplied by the coefficient \(\beta_1\). \(x\) is dummy-coded, where indicates that the observed \(y_i\) is an effect size estimate of compared to control, while indicates that \(y_i\) is an effect size estimate of the compared to control.
Here are the priors for this comparison:
\[ \begin{align*} \beta_0 & \sim \operatorname{Normal}(0, 0.75^2) \tag{Priors}\\ \tau & \sim \operatorname{Log-Normal}(-1.975, 0.67^2)\\ \beta_{1[Skeptical]} & \sim \operatorname{Normal}(0, 0.354^2)\\ \beta_{1[Sarilumab]} & \sim \operatorname{Normal}(-0.049, 0.193^2)\\ \beta_{1[Tocilizumab]} & \sim \operatorname{Normal}(0.049, 0.193^2)\\ \beta_{1[Vague]} & \sim \operatorname{Normal}(0, 4^2) \end{align*} \]
## Means
# Figure S3 in https://doi.org/10.1101/2021.06.18.21259133;
# version posted June 22, 2021
reciprocal_remap_cap = 1/1.05
mean_log_sari = log(reciprocal_remap_cap)
mean_skeptical = log(1)
## SDs
# Assuming a mean = log(1), what is the SD that yields a distribution
# with 0.025 probability density below log(0.5)?
sari_sd_skeptical = (mean_skeptical - log(0.5))/qnorm(1-0.025)
# Assuming a mean = log(0.952), what is the SD that yields a distribution
# with 0.4 probability density above log(1)?
sari_sd_optimistic_sari = (log(1) - mean_log_sari)/qnorm(1-0.4)
# Assuming a mean = log(1/0.952) = -log(0.952), what is the SD that yields a
# distribution with 0.4 probability density below log(1)?
sari_sd_optimistic_toci = (-mean_log_sari - log(1))/qnorm(1-0.4)
sari_sd_vague = 4
# General model characteristics
## Formula
mf =
formula(yi | se(sei) ~ 0 + Intercept + treatment + (1 | study))
## Priors for beta_0 and tau
# tau:
# logmean
# -1.975
# logsd
# 0.67
priors =
brms::prior(normal(0, 0.75), class = "b", coef = "Intercept") +
brms::prior(lognormal(-1.975, 0.67), class = "sd", group = "study")
Fit models
## Skeptical model
# sari_sd_skeptical
# 0.353653
m_sari_skeptical =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "baricitinib"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(0, 0.353653), class = "b", coef = "treatmentsarilumab"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_sari_skeptical.Rds"),
file_refit = "on_change"
)
## "Optimistic for Sarilumab" model
# mean_log_sari
# -0.04879016
# sari_sd_optimistic_sari
# 0.1925823
m_sari_optimistic_sari =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "baricitinib"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(-0.04879016, 0.1925823), class = "b", coef = "treatmentsarilumab"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_sari_optimistic_sari.Rds"),
file_refit = "on_change"
)
## "Optimistic for Tocilizumab" model
# -mean_log_sari
# 0.04879016
# sari_sd_optimistic_toci
# 0.1925823
m_sari_optimistic_toci =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "baricitinib"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(0.04879016, 0.1925823), class = "b", coef = "treatmentsarilumab"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_sari_optimistic_toci.Rds"),
file_refit = "on_change"
)
## Vague model
m_sari_vague =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "baricitinib"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(0, 4), class = "b", coef = "treatmentsarilumab"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_sari_vague.Rds"),
file_refit = "on_change"
)
Here are the priors for this comparison:
\[ \begin{align*} \beta_0 & \sim \operatorname{Normal}(0, 0.75^2) \tag{Priors}\\ \tau & \sim \operatorname{Log-Normal}(-1.975, 0.67^2)\\ \beta_{1[Skeptical]} & \sim \operatorname{Normal}(0, 0.354^2)\\ \beta_{1[Baricitinib]} & \sim \operatorname{Normal}(-0.105, 0.416^2)\\ \beta_{1[Tocilizumab]} & \sim \operatorname{Normal}(0.105, 0.416^2)\\ \beta_{1[Vague]} & \sim \operatorname{Normal}(0, 4^2) \end{align*} \]
## Means
mean_log_bari = log(0.9)
mean_skeptical = log(1)
## SDs
# Assuming a mean = log(1), what is the SD that yields a distribution
# with 0.025 probability density below log(0.5)?
bari_sd_skeptical = (mean_skeptical - log(0.5))/qnorm(1-0.025)
# Assuming a mean = log(0.9), what is the SD that yields a distribution
# with 0.4 probability density above log(1)?
bari_sd_optimistic_bari = (log(1) - mean_log_bari)/qnorm(1-0.4)
# Assuming a mean = log(1/0.9) = -log(0.9), what is the SD that yields a
# distribution with 0.4 probability density below log(1)?
bari_sd_optimistic_toci = (-mean_log_bari - log(1))/qnorm(1-0.4)
bari_sd_vague = 4
## Skeptical model
# bari_sd_skeptical
# 0.353653
m_bari_skeptical =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(0, 0.353653), class = "b", coef = "treatmentbaricitinib"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_bari_skeptical.Rds"),
file_refit = "on_change"
)
## "Optimistic for Baricitinib" model
# mean_log_bari
# -0.1053605
# bari_sd_optimistic_bari
# 0.4158742
m_bari_optimistic_bari =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(-0.1053605, 0.4158742), class = "b", coef = "treatmentbaricitinib"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_bari_optimistic_bari.Rds"),
file_refit = "on_change"
)
## "Optimistic for Tocilizumab" model
# -mean_log_bari
# 0.1053605
# bari_sd_optimistic_toci
# 0.4158742
m_bari_optimistic_toci =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(0.1053605, 0.4158742), class = "b", coef = "treatmentbaricitinib"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_bari_optimistic_toci.Rds"),
file_refit = "on_change"
)
## Vague model
m_bari_vague =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(0, 4), class = "b", coef = "treatmentbaricitinib"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_bari_vague.Rds"),
file_refit = "on_change"
)
## Skeptical model
# bari_sd_skeptical
# 0.353653
m_bari_skeptical_sens =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari (No Toci)"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(0, 0.353653), class = "b", coef = "treatmentbaricitinib"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_bari_skeptical_sens.Rds"),
file_refit = "on_change"
)
## "Optimistic for Baricitinib" model
# mean_log_bari
# -0.1053605
# bari_sd_optimistic_bari
# 0.4158742
m_bari_optimistic_bari_sens =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari (No Toci)"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(-0.1053605, 0.4158742), class = "b", coef = "treatmentbaricitinib"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_bari_optimistic_bari_sens.Rds"),
file_refit = "on_change"
)
## "Optimistic for Tocilizumab" model
# -mean_log_bari
# 0.1053605
# bari_sd_optimistic_toci
# 0.4158742
m_bari_optimistic_toci_sens =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari (No Toci)"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(0.1053605, 0.4158742), class = "b", coef = "treatmentbaricitinib"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_bari_optimistic_toci_sens.Rds"),
file_refit = "on_change"
)
## Vague model
m_bari_vague_sens =
brms::brm(
data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari (No Toci)"),
family = gaussian,
formula = mf,
prior = priors +
brms::prior(normal(0, 4), class = "b", coef = "treatmentbaricitinib"),
sample_prior = TRUE,
control = list(adapt_delta = .96),
backend = "cmdstanr", # faster
warmup = 2000,
iter = 4000,
seed = 123,
cores = parallel::detectCores(),
chains = 4,
file = here::here("output", "fits", "mr_bari_vague_sens.Rds"),
file_refit = "on_change"
)
# Load function
source(here::here("functions", "diag_plot.R"))
Tocilizumab model
diag_plot(model = ma_toci,
pars_list = c("b_Intercept", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(ma_toci,
type = "dens_overlay",
ndraws = 50)
brms::pp_check(ma_toci,
type = "ecdf_overlay",
ndraws = 50)
Sarilumab model
diag_plot(model = ma_sari,
pars_list = c("b_Intercept", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(ma_sari,
type = "dens_overlay",
ndraws = 50)
brms::pp_check(ma_sari,
type = "ecdf_overlay",
ndraws = 50)
Baricitinib models
diag_plot(model = ma_bari,
pars_list = c("b_Intercept", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(ma_bari,
type = "dens_overlay",
ndraws = 50) + coord_cartesian(y = c(0, 5))
brms::pp_check(ma_bari,
type = "ecdf_overlay",
ndraws = 50)
Baricitinib models
diag_plot(model = ma_bari_sens,
pars_list = c("b_Intercept", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(ma_bari_sens,
type = "dens_overlay",
ndraws = 50) + coord_cartesian(y = c(0, 5))
brms::pp_check(ma_bari_sens,
type = "ecdf_overlay",
ndraws = 50)
“Skeptical” model
diag_plot(model = m_sari_skeptical,
pars_list = c("b_Intercept", "b_treatmentsarilumab", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_sari_skeptical,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50)
brms::pp_check(m_sari_skeptical,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
“Optimistic for Sarilumab” model
diag_plot(model = m_sari_optimistic_sari,
pars_list = c("b_Intercept", "b_treatmentsarilumab", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_sari_optimistic_sari,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50)
brms::pp_check(m_sari_optimistic_sari,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
“Optimistic for Tocilizumab” model
diag_plot(model = m_sari_optimistic_toci,
pars_list = c("b_Intercept", "b_treatmentsarilumab", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_sari_optimistic_toci,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50)
brms::pp_check(m_sari_optimistic_toci,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
“Vague” model
diag_plot(model = m_sari_vague,
pars_list = c("b_Intercept", "b_treatmentsarilumab", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_sari_vague,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50)
brms::pp_check(m_sari_vague,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
“Skeptical” models
diag_plot(model = m_bari_skeptical,
pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_bari_skeptical,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50) + coord_cartesian(y = c(0, 4))
brms::pp_check(m_bari_skeptical,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
diag_plot(model = m_bari_skeptical_sens,
pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_bari_skeptical_sens,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50) + coord_cartesian(y = c(0, 4))
brms::pp_check(m_bari_skeptical_sens,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
“Optimistic for Baricitinib” models
diag_plot(model = m_bari_optimistic_bari,
pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_bari_optimistic_bari,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50) + coord_cartesian(y = c(0, 4))
brms::pp_check(m_bari_optimistic_bari,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
diag_plot(model = m_bari_optimistic_bari_sens,
pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_bari_optimistic_bari_sens,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50) + coord_cartesian(y = c(0, 4))
brms::pp_check(m_bari_optimistic_bari_sens,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
“Optimistic for Tocilizumab” models
diag_plot(model = m_bari_optimistic_toci,
pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_bari_optimistic_toci,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50) + coord_cartesian(y = c(0, 4))
brms::pp_check(m_bari_optimistic_toci,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
“Optimistic for Tocilizumab” model
diag_plot(model = m_bari_optimistic_toci_sens,
pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_bari_optimistic_toci_sens,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50) + coord_cartesian(y = c(0, 4))
brms::pp_check(m_bari_optimistic_toci_sens,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
“Vague” models
diag_plot(model = m_bari_vague,
pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_bari_vague,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50) + coord_cartesian(y = c(0, 4))
brms::pp_check(m_bari_vague,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)
diag_plot(model = m_bari_vague_sens,
pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
ncol_trace = 4)
brms::pp_check(m_bari_vague_sens,
type = "dens_overlay_grouped",
group = "treatment",
ndraws = 50) + coord_cartesian(y = c(0, 4))
brms::pp_check(m_bari_vague_sens,
type = "ecdf_overlay_grouped",
group = "treatment",
ndraws = 50)